Ab initio calculation of excitonic effects in the optical spectra of semiconductors 
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An ab initio approach to the calculation of excitonic effects in the optical absorption spectra of 
semiconductors and insulators is formulated. It starts from a quasiparticle bandstructure calculation 
and is based on the relevant Bethe-Salpeter equation. An application to bulk silicon shows a 
substantial improvement with respect to previous calculations in the description of the experimental 
spectrum, for both peak positions and lineshape. 
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Recent advances in ab initio calculations, mostly Den- 
sity Functional Theory - Local Density Approximation 
(DFT-LDA) applications, allow to determine the ground 
state properties and the Kohn-Sham (KS) electronic 
structure Jl| for even complicated systems. In order to 
treat excited states, realistic quasiparticle (QP) energies 
are then in general obtained by applying self-energy cor- 
rections to the KS energies, usually evaluated in the GW 
approximation Excellent agreement of the resulting 
bandstructure with experimental data has been found 
for a wide range of materials ||,^|]. However, spectro- 
scopic properties involving two-particle excitations are 
often only poorly described at this one-particle level. The 
main example is absorption spectroscopy, where a simul- 
taneously created electron-hole pair interacts more or 
less strongly. As a consequence, in addition to bound 
exciton states which occur within the gap, the spectral 
lineshape above the continuous-absorption edge is dis- 
torted. 

The reported qualitative agreement with experiment 
of many computed KS-LDA absorption spectra, obtained 
from one-electron transitions between KS states |5j , is in- 
deed due to a partial cancellation between two principal 
errors: namely, the compensation of the large KS-LDA 
underestimation of the valence-conduction bandgap, 
with an overestimation of the absorption onset induced 
by calculating the dielectric function entirely within the 
one-particle picture. The situation often worsens when 
only the first error is corrected by replacing the KS eigen- 
values with the realistic QP energies |6[0]. On the other 
hand, going beyond the one-particle picture through in- 
clusion of local field and/or exchange-correlation effects 
within DFT-LDA in the calculation of the absorption 
spectrum does generally not remove the observed dis- 
crepancy ||. In fact, most of the residual error stems 
from the neglect of the electron-hole interaction. 
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Up to now, excitonic effects have been rarely calculated 
from first principles. Some information about energetic 
changes can be extracted from an LDA-based A-SCF ap- 
proach M. Large excitonic effects on the spectral prop- 
erties have been calculated ab initio in the case of a small 
sodium cluster || . This approach has consequently been 
generalized with the calculation of the absorption onset 
for an infinite system, the Li 2 crystal JhJ. The calcu- 
lation of the entire optical spectrum of a solid, finally, 
still remains a major challenge Quantitatively cor- 
rect theoretical absorption spectra are indeed needed as 
a reference for the interpretation and prediction of ex- 
perimental results. 

A paradigmatic case is bulk silicon, which is represen- 
tative for the group IV, III-V, and II- VI semiconductors. 
These materials show qualitatively similar optical spec- 
tra, with two major structures at 3-5 eV. The first peak 
(Ei) has been interpreted as a Mi type critical point 
transition, and the second peak (E 2 ) as a M 2 type one 
p3[ . Theoretical work based on the one-electron approxi- 
mation, ranging from early empirical pseudopotential ap- 
proaches 1 14 1 to ab initio DFT-LDA work all yielded 
the same qualitative result, i.e. an underestimation of 
the Ei peak by as much as 50%, reducing it to a weak 
shoulder of the generally overestimated E 2 peak. In order 
to go beyond, Louie et al. included local field effects 
in the calculation of the dielectric matrix. The result- 
ing spectrum is significantly improved at higher energies 
(above 15 eV), but not in the region of interest around 4 
eV. 

Several authors suggested that strong contributions 
to the Ei peak could arise from saddle point excitons 
pG| |i"8[ . Excitonic effects allowed to explain the mea- 
sured temperature and pressure dependence of the line- 
shape and the symmetry in wavelength modulation re- 
flectance spectra |l7j . Until now, the most sophisticated 
calculation of excitonic effects on the spectral lineshape 
of silicon was done by Hanke and Sham |l8| . They per- 
formed a semi-empirical LCAO calculation, including lo- 
cal field effects and the screened electron-hole attraction. 
As in Ref. ifTsf , local field effects alone were shown to 
transfer oscillator strength to higher energies and hence 
to increase the discrepancy with experiment at lower en- 
ergies. On the contrary, the electron-hole interaction 
shifted the position of the Ei peak to lower energies, and 
almost doubled its intensity, while the oscillator strength 
of the higher energy peaks was decreased. The over- 
all agreement with experiment was hence improved, and 
clear evidence was given for the importance of excitonic 
effects. However, the final intensity ratio between the 
Ei and E 2 peaks was reversed, in disagreement with the 
experimental spectrum. As pointed out by Wang et al. 
p3[ , the reliability of semi-empirical approaches is lim- 
ited. For instance, there are important differences, al- 
ready at the one-electron level, between the spectra of 
Refs. H and @. 

In the ab initio framework, on the other hand, the pre- 
cision achievable for the computation of electronic spec- 
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tra is in general still poor when compared with the qual- 
ity of calculated ground state properties. This work is 
aimed to shrink this gap, showing how a significant im- 
provement of the ab initio calculation of absorption spec- 
tra can be obtained. 

The absorption spectrum is given by the imaginary 
part of the macroscopic dielectric function cm 

e M (u) = 1 - lim v(q)XG=o,G'=o(q;w), W 

q^O 

where x(r, r'; u) = -iS(r, r, r', r'; u>). 5(1, T; 2, 2') is the 
part of the two-particle Green's function which excludes 
the disconnected term -G(l, l')G(2, 2'), and G(l, 1') is 
the one-particle Green's function p8| . The notation 
(1,2) stands for two pairs of space and time coordinates, 
(r 1 ,t 1 ;r 2 ,t 2 ). 

Following Ref. fl8|| , we start from the Bethe-Salpeter 
equation for 5, 

5(1, 1'; 2, 2') = 5 (1, 1'; 2, 2') + 5 (1, 1'; 3, 3')S(3, 3'; 4, 4')5(4, 4'; 2, 2'). 

(2) 

Repeated arguments are integrated over. The term 
5 (1,1';2,2') = G(1',2')G(2,1) yields the polarization 
function of independent quasiparticles \q, from which 
the standard RPA €m is obtained. The kernel S contains 
two contributions: 

5(1, 1', 2, 2') - -iS(l, 1>(2, 2>(1, 2) + iS(l, 2)6(1', 2')W(1, l'). 

(3) 

Considering the first term in the calculation of S is 
equivalent to the inclusion of local field effects in the ma- 
trix inversion of a standard RPA calculation. In order 
to obtain the macroscopic dielectric constant, the bare 
Coulomb interaction v contained in this term must, how- 
ever, be used without the long range term of vanishing 
wave vector [|l9|. When spin is not explicitly treated, 
v gets a factor of two for singlet excitons. In the sec- 
ond term, W is the screened Coulomb attraction between 
electron and hole. It is obtained as a functional deriva- 
tive of the self-energy in the GW approximation, neglect- 
ing a term G^j. This latter term contains information 
about the change in screening due to the excitation, and 
is expected to be small [^(J . We limit ourselves to static 
screening, since dynamical effects in the electron-hole 
screening and in the one particle Green's function tend 
to cancel each other M] , which suggests to neglect both 
of them. 

In order to solve Eq. (2), we have to invert a 4-point 
function. In Ref. |Q this has been possible due to the 
use of a very limited basis set. In an ab initio plane wave 
calculation, such a procedure is clearly prohibitive, when 
plane waves are chosen as straightforward basis functions. 
Instead, the physical picture of interacting electron-hole 
pairs suggests to use a basis of LDA Bloch functions, 
ip n (r), expecting that only a limited number of electron- 
hole pairs will contribute to each excitation. 
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In this basis, ^"^"^(^ru) _ 5 nun3 S n2 ^ ni (f n2 - 
f ni )/(E n2 — E ni — uj) and, after solving for S, in the 
case of static screening, equation (2) can be written as 

^ r (ni,n2),(n 3 ,Ti4) = (Hexc ~ ^ w )( ni , r , 2 ),(„ 3 ^ 4 )(/n4 ~ /n 3 )> 

(4) 

with 

rr(ni,n 2 ),(n 3 ,n 4 ) _ ( J? _R U A — iff —f Ix 

1 dn J dr[J dr 2 J dr' 2 ^ (n) ^(rj) H( ri , ri, r 2 , r 2 ) ^ 3 (r 2 ) ^ 4 (r 2 ). (5) 

J is the identity operator. The energies -E„ are the QP 
levels. Together with the above form of xo this is con- 
sistent with the use of LDA wavefunctions and updated 
energy denominators in the Green's function used to con- 
struct the self-energy in the GW calculation. The f n are 
Fermi-Dirac occupation numbers. We avoid to invert the 
matrix (H exc — I uj) for each absorption frequency u> by 
applying the identity 

. ^ |A > M7\, < A'| 
(H exc -I^= ^ (E » u) (6) 

which holds for a system of eigenvectors and eigenvalues 
of a general, non-hermitian matrix defined by 

H exc \X>=E x \X> . (7) 

M\ \i is the overlap matrix of the (in general non- 
orthogonal) eigenstates of H exc . 

Equation (7) is the effective two-particle Schrodinger 
equation which we solve by diagonalization. The explicit 
knowledge of the coupling of the various two-particle 
channels, given by the coefficients A^ 1 '™ 2 ^ of the state 
| A > in our LDA basis, allows to identify the character 
of each transition. (This analysis would be much more 
cumbersome if a matrix inversion instead of the spectral 
representation was chosen, as in Ref. |18|.) 

The macroscopic dielectric function in Eq. (1) is ob- 
tained as 

e M (u) = 1 - lim «(q) £ M^{, ]T < rn^'fa > x 

A, A' Til .712 



E 



<n 4 |e + ^ r |n 3 >^" 3 '" 4) ( {" 4 fn *} . (8) 

{bj\ - LU) 



In practice, the KS eigenvalues and eigenfunctions from 
a DFT-LDA calculation serve as input to the evaluation 
of the RPA screened Coulomb interaction W and the 
GW self-energy E. The KS eigenfunctions, together with 
the QP energies and W, are then used in the exciton 
calculation. Here each pair of indices (ni,n 2 ) stands for 
a pair of bands and one Bloch vector k in the Brillouin 
zone (BZ), since we are interested in direct transitions 
only. 
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In principle, all combinations of bands should be con- 
sidered. It can, however, be shown exactly that only pairs 
containing one filled and one empty LDA state contribute 
to (8). Still, the portion of the matrix H exc to be consid- 
ered is in general non-hermitian, being of the form p2| 



H 



Jj(v 1 ,c 1 ),(v 2 ,c 2 ) H^ Vl ' Cl )>( C2 '' U2 ) 
_JJ*(v 1 .c 1 ),(c 2 ,v 2 ) _JJ*(v 1 ,c 1 ),(v 2 ,c 2 ) 



The off-diagonal coupling matrices do not contain the 
QP transition energies, but only the interaction elements, 
which arc much smaller in the case of silicon. Hence, we 
neglect the latter and separate the Hamiltonian into two 
block-diagonal parts: the resonant contributions, which 
are active for positive frequencies, and the antiresonant 
ones, only contributing to negative frequencies. The ma- 
trix of the resonant part by its own is hermitian, and we 
therefore obtain the simpler formula 



emM = 1 + lim v(q) Y] 

q^O ^ — ' 



S,, c;k <^k|e-^|c,k>4^ k) | 2 
{E x - u) 

(9) 



(7) and (9) constitute a set of equations which has been 
frequently used in the non-a& initio framework pp| , p4[ . 
Here, it appears as a particular approximation to the 
more general formula (8), with well-defined ab initio in- 
gredients which are consistent with the GW approach. 

We evaluate expression (9) for bulk silicon. The DFT- 
LDA calculation is performed using norm-conserving 
pseudopotentials |2q| , an energy cutoff of 15 Ry, and 256 
special k points in the BZ (^6|. Next, GW corrections 
to the KS band structure are obtained following the ap- 
proach of j27|. The quite smooth GW corrections arc 
interpolated for the denser k point mesh needed for the 
absorption spectrum. We evaluate equations (7) and (9) 
using different sets containing up to 2048 k points in the 
BZ. In order to handle such large matrices, the symme- 
try properties of the crystal are exploited. One has to be 
very careful in doing so, since the spectrum turns out to 
be extremely sensitive to any inconsistency in the phases 
which may appear when wavefunctions are rotated, no- 
tably for degenerate bands. A safe way to proceed is 
to make only partial use of symmetry, considering only 
those operations which form an abelian subgroup of the 
point group, and which altogether allow to reconstruct 
the whole zone from a corresponding reduced part. In 
the case of silicon, we found it convenient to use the 
180° rotations around the x and the y axis, respectively. 
These two operations T allow us to break the equation 
H A k = EA k (band indices have been suppressed, 
and repeated indices are summed over) into four equa- 
tions to be used for points k^ in the reduced zone only. 
These equations are of the form h^k^a 1 = Ea ki , where 
/i k . k / . = H kik i ±H kiTk i . The A k are then reconstructed 
from the reduced eigenvectors a ki . Moreover, we apply 
time reversal and hermiticity in order to accelerate the 
calculation of the matrix elements. 
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A set with 864 k points in the full BZ is used to 
check the various ingredients of our calculation, in par- 
ticular the number of bands and the importance of the 
off-diagonal elements of the inverse dielectric matrix in 
the evaluation of W. In the inset of Fig. 1, the continu- 
ous line shows the results of a calculation with 4 valence 
and 4 conduction bands, and the full e _1 . In the re- 
gion of interest (below 4.5 eV), a 6 bands calculation (4 
valence + 2 conduction, dotted line) appears to be suf- 
ficient. Neglecting the off-diagonal elements of £qq, (q) 
yields an indistinguishable curve. We then use 6 bands 
and the diagonal (T 1 to compute the spectrum with 2048 
k points in the full BZ. In the main part of Fig. 1 the ex- 
perimental spectrum (dotted line) j28| is compared to: i) 
an RPA calculation |]2j| taking only into account the QP 
shifts, but not the excitonic or local field effects (short- 
dashed curve): the result is, as generally observed, in 
great discrepancy with experiment; ii) a calculation in- 
cluding local field effects (i.e. using equation (3) with W 
set to 0, long-dashed curve): the agreement is worsened, 
since the oscillator strength is slightly shifted to higher 
energies and both the Ex and E 2 peaks are lowered, thus 
confirming previous findings in the literature |l5| , |l8| ] ; iii) 
finally, the full calculation including the electron-hole 
attraction (continuous curve): absolute intensities now 
agree well with experiment. The remaining slight overes- 
timate is of the order of what has been predicted by Ref. 
|2lf to be the contribution of dynamical effects. More 
important, the peak positions and the relative intensity 
of the main structures are both in good agreement with 
experiment. Also the structure at 3.8 eV, even though 
slightly overestimated due to a finite k point sampling, 
has been repeatedly confirmed in theoretical and experi- 
mental work ||(J . 

In conclusion, we have shown how excitonic effects can 
be included in an ab initio calculation of optical absorp- 
tion spectra of semiconductors. At the example of bulk 
silicon, we have demonstrated that good agreement with 
experiment can be obtained for a case where the inclusion 
of self-energy and local field effects alone still gives rise to 
a rather poor theoretical spectrum. In this context, bulk 
silicon is not particularly easy to handle, since the bottle- 
neck of the calculation is the number of k points (high in 
silicon, due to large dispersion) and not the energy cutoff. 
Even though, the computational effort, mostly steming 
from diagonalization, is reasonable and demands only a 
few hours on a Cray C98. The present work opens hence 
the way to first-principles calculations of optical absorp- 
tion spectra with a precision comparable to that typically 
achieved in ground state calculations. 
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FIG. 1. Absorption spectra of Si. Inset: calculation ac- 
cording to equation (9) with 864 k points in the BZ, using 
8 bands (continuous curve) or only 6 bands (dotted curve). 
Main part: Calculation according to equation (9) with 2048 k 
points in the BZ, 6 bands and the diagonal approximation to 
e _1 : with both electron-hole attraction and local field effects 
in the Hamiltonian (continuous curve), inclusion of local field 
effects alone (long-dashed curve) and RPA with QP shifts 
only (short-dashed curve). Experimental curve (dots) [28]. 
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